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INTRODUCTION 

This  program  builds  off  of  our  extensive  experience  in  using  electrical  properties  of  prostate  to  distinguish 
malignant  from  benign  tissues  [1-5]  and  specifically  stems  from  exciting  new  data  published  in  The  Prostate  [6] 
in  which  we  demonstrated  significant  electrical  property  differences  between  high-  and  low-grade  prostate 
cancer.  These  electrical  properties  are  influenced  by  a  tissue’s  intra-  and  extra-cellular  composition, 
morphology,  and  cellular  constituency,  and  we  have  hypothesize  that  it  is  possible  to  use  these  properties  to 
discriminate  between  normal,  low-grade,  and  high-grade  malignant  formations  in  a  clinical  setting.  While 
measuring  these  properties  by  direct  contact  with  the  tissues  is  possible  in  invasive  experiments,  it  is  desirable 
to  develop  methods  to  do  so  in  a  non-invasive  fashion.  To  date  our  group,  and  other  groups  around  the  world, 
have  investigated  Electrical  Impedance  Tomography  and  Microwave  Imaging,  two  techniques  which  are  limited 
in  resolution  by  the  underlying  physics.  The  primary  objective  of  this  current  program  is  to  develop  a  high- 
resolution  MR-based  approach  to  imaging  the  electrical  properties  of  prostate  with  the  intent  of  producing  a 
system  potentially  able  to  image  cancer  grade.  This  is  possible  by  leveraging  ultra-novel  developments  in  MRI, 
and  our  extensive  experience  in  developing  technologies  to  gauge  and  assess  the  utility  of  electrical  properties 
for  prostate  cancer  detection  [1-6]  and  assessment  and  in  developing  computer  algorithms  to  transform 
electromagnetic  data  into  electrical  property  images  of  the  prostate  [7-14].  Specifically,  we  are  attempting  to 
use  the  maps  of  MRI  RF  field  data  acquired  with  safe  and  fast  sequences  to  create  high-resolution  electrical 
property  images  of  the  prostate.  We  are  developing  this  novel  technology,  evaluating  it  in  an  ex  vivo  setting, 
and  finally  assessing  the  feasibility  of  employing  this  imaging  modality  in  a  routine  clinical  cohort  of  patients 
with  the  intent  of  having  a  significant  and  immediate  impact  on  clinical  practice.  By  developing  this  high- 
resolution  electrical  property  imaging  modality  we  expect  to  produce  highly  sensitive  and  specific  images  of 
cancer  grade  within  the  prostate  and  ultimately  better  guide  clinicians  in  distinguishing  aggressive  from 
indolent  disease. 

Much  of  the  third  year  of  this  program  has  focused  on  ex  vivo  prostate  imaging,  ex  vivo  data  analysis, 
continued  developing  of  MR-EPT  conductivity  reconstructions,  and  initiation  of  in  vivo  data  collections.  We  are 
currently  in  a  No  Cost  Extension  period  (additional  6  months)  in  which  we  will  focus  on  completing  the  in  vivo 
data  collection  and  statistical  analysis  of  our  data.  In  addition  we  will  be  preparing  publications  and  proposals 
focused  more  heavily  on  clinical  data  acquisition  and  evaluation. 

BODY 

The  following  research  summary  is  presented  in  terms  of  the  approved  Statement  of  Work,  with  each  task 
being  discussed  separately.  When  appropriate,  detailed  discussion  is  referenced  to  manuscripts  published, 
submitted,  or  in  preparation  which  are  provided  in  the  Appendix.  Information  provided  in  our  previous  annual 
report  is  omitted  here  and  instead  we  note  the  tasks  that  are  completed  and  reference  our  2014/2015  annual 
reports.  Note  that  future  task  and  objectives  to  be  completed  are  marked  as  TBC. 

SPECIFIC  AIM  1:  TO  DEVELOP  MR-EPT  FOR  PROSTATE  IMAGING 
Major  Task  1 :  Develop  computation  toolbox  for  MR-EPT 

a)  Build  Matlab-based  toolbox  for  computing  electrical  property  images 

Completed  and  reported  on  in  our  2014/2015  Annual  Reports.  On-going  investigation  continue  in  order  to 
further  optimize  and  validate  the  approach  developed. 

We  have  fully  implemented  our  proposed  method  as  a  MATLAB  toolbox  for  MR-EPT  image  reconstruction 
as  described  in  our  2014/2015  annual  reports.  Over  the  past  year  we  made  two  significant  advancements. 
The  first  is  a  method  for  smoothing  the  boundaries  between  sub-domains  within  our  images  and  the 
second  is  a  method  for  reconstructing  permittivity  of  an  object.  In  terms  of  smoothing  out  the  sub-domains, 
we  have  developed  an  overlapping  method  by  which  we  reconstruct  conductivity  over  multiple  overlapping 
sub-domains  of  different  sizes.  By  properly  accounting  for  the  overlapping  pixels  we  can  average  these 
sub-domains  together  and  smooth  out  the  boundaries  between  the  sub-domains.  This  creates  significantly 
more  appealing  images  that  do  not  have  sharp  transitions  between  the  different  sub-domains  (as  were 
present  in  our  previous  implementations).  Figure  1  depicts  the  process  and  Figure  2  demonstrates  the 
significant  image  quality  improvements  achievable  with  this  optimized  algorithm. 


Figure  1.  Illustration  of  the  tiling  subdomains,  data  used  in  reconstructing  (red  square),  and  a  2nd  derivative  smoothing  constrain. 
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Figure  2:  Several  examples  of  multiple  overlays  used  to  smooth  the  internal  boundary  artifacts.  Note  that  the  image  on  the  left  was  our 
initial  implementation  and  that  the  internal  artifacts  are  significantly  reduced  as  we  average  over  multiple  overlays 


In  terms  of  permittivity  imaging,  we  have  developed  a  forward  model,  the  Jacobian,  and  an  iterative 
approach  to  estimating  the  permittivity  distribution.  We  have  evaluated  this  in  simulation  and  plan  to 
evaluate  this  in  phantom  models  during  the  next  quarter.  Figure  3  shows  an  example  image  of  permittivity 
reconstruction  and  a  PowerPoint  presentation  describing  this  development  is  appended  to  this  quarterly 
report. 
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Figure  3:  Example  simulation  of  permittivity  reconstruction.  Left  image  is  true  permittivity  distribution.  Central  image  represents  the  B1 
amplitude  map  computed  with  our  forward  solver.  Right  image  represents  the  reconstructed  permittivity.  Note  that  the  reconstructed 
permittivity  closely  matches  the  true  permittivity  distribution  (left  image),  with  the  exception  of  a  few  boundary  pixels  which  are 
artificially  elevated  due  to  how  the  boundary  conditions  are  implemented  to  solve  the  forward  problem. 


b)  Build  Matlab-based  toolbox  for  specific  MR-based  field  of  views 

Completed  and  reported  on  in  our  2014/2015  Annual  Reports.  On-going  investigation/developments 
continue  in  order  to  further  optimize  and  validate  this  toolbox. 

We  have  developed  Matlab-based  functions  to  read  in  arbitrary  MRI  DICOM  and  .PAR/.REC  (Phillips 
format)  images  and  display  them  for  evaluation.  Additional  developments  have  included  the  ability  to  create 
multi-slice  images  of  the  MR  data  for  use  in  comparing  the  different  MR  variants  we  are  exploring.  This 
toolbox  was  briefly  described  in  our  2014/2015  annual  reports. 

Task  1  Milestones: 

1 .  Functional  toolbox  for  producing  MR-EPT  images  -  Completed 


Major  Task  2:  Optimize  MR-EPT  and  multi-parametric  MR  imaging  through  phantom  imaging 

a)  Optimize  MR-sequence  for  MR-EPT 

Completed  and  reported  on  in  our  2014/2015  Annual  Reports. 

b)  Perform  initial  tank-based  phantom  imaging  studies 

Completed  and  reported  on  in  our  2014/2015  Annual  Reports.  On-going  investigation  continue  in  order  to 
further  optimize  and  validate  the  approach  developed. 

A  number  of  additional  tank-based  phantom  studies  were  conducted  during  this  past  year  to  further 
demonstrate  that  our  MR-EPT  algorithms  are  able  to  accurately  estimate  the  internal  conductivity  of  a 
volume.  One  such  experiment  specifically  focused  on...  validating  that  our  approach  is  able  to  accurately 
reconstruct  curved  surfaces  (some  of  our  previous  experiments  explored  linearly  varying  conductivity 
changes).  Figure  2  displays  a  curvilinear  gelatin  phantom  conducted  along  with  the  MR  magnitude  image 
produced.  The  phase  data  recorded  from  this  experiment  was  reconstructed  using  both  QR  and  TV 
approaches  to  solve  the  inverse  problem  (Figure  3).  Both  approaches  accurately  reconstruct  the  curved 
surface  suggesting  that  these  techniques  extend  beyond  reconstructing  linearly  varying  conductivity 
distributions.  In  addition,  we  divide  the  phase  images  into  a  number  of  sub-domains  to  reduce  the 
computational  time  associated  with  image  reconstruction.  We  explored  this  and  the  time  associated  with 
using  different  numbers  of  sub-domains  (Figure  3  and  Table  1).  The  experimental  configuration  and 
analysis  are  further  described  in  Appendix  1. 

c)  Perform  anatomically  accurate  phantom  imaging  studies 
Completed  and  reported  on  in  our  2014/2015  Annual  Reports. 

Task  2  Milestones: 

Note  that  while  the  below  tasks  are  complete,  we  expect  to  continue  conducting  phantom  experiments  over  the 
course  of  the  next  year  to  continue  to  improve  our  image  reconstruction  algorithms  and  better  understand  any 
image  artifacts  that  may  appear  in  our  clinical  data  acquisition. 

1 .  Validated  MR-EPT  algorithms  -  Completed 

2.  Fully  functional  protocol  for  obtaining  MR-based  images  in  a  single  serially  acquired  imaging  session  - 

Completed 

3.  1  peer-reviewed  publication  submitted  -  Completed,  see  Appendix 

Major  Task  3:  Submit  documents  for  IRB  and  MRMC  HRPO  approval 

a)  Draft  and  submit  IRB  protocol  revisions  and  new  protocol  submission 
Ex  vivo  Protocol:  Completed  during  last  annual  reporting  period. 

In  Vivo  Protocol:  Completed  during  last  annual  reporting  period 

b)  Draft  and  submit  documentation  for  MRMC  HRPO  approval 

Ex  vivo  Protocol:  Completed  during  last  annual  reporting  period 

In  Vivo  Protocol:  Completed  during  last  annual  reporting  period 

Task  3  Milestones: 

1.  Obtain  IRB  and  MRMC  HRPO  approval  for  ex  vivo  and  in  vivo  cohorts  -  ex  vivo  completed,  in  vivo 
completed 

SPECIFIC  AIM  2:  TO  EVALUATE  MR-EPT  IN  AN  EX  VIVO  COHORT  OF  PROSTATES 
Major  Task  1:  Optimize  ex  vivo  MR-EPT  and  multi-parametric  MR  imaging 

a)  Record  MR-EPT  and  multi-parametric  MR  sequences  of  ex  vivo  prostates 
This  task  has  been  complete  and  described  in  our  2014/2015  Annual  Reports. 

b)  Optimize  MR-EPT  sequences  and  algorithms  based  on  findings  in  this  initial  cohort 
This  task  has  been  complete  and  described  in  our  2014/2015  Annual  Reports. 


Task  1  Milestones: 

1.  Validation  that  our  MR-EPT  and  multi-parametric  MR  protocol  is  initially  optimized,  robust,  and 
repeatable  -  Completed 

Major  Task  2:  Evaluate  ex  vivo  MR-EPT  and  multi-parametric  MR  imaging 

a)  Record  multi-parametric  MR  sequences  of  ex  vivo  prostates 

Over  the  past  year  we  have  actively  recruited  patients  to  participate  in  our  ex  vivo  study.  We  have  been 
averaging  ~  2  cases  per  month.  To  date  we  have  imaged  50  ex  vivo  prostates  (reaching  our  target 
accrual).  For  all  cases  we  have  been  recording  B1  phase  and  magnitude  images  (for  MR-EPT  image 
reconstruction),  T1  spin  echo,  T1  turbo  spin  echo,  and  T2  turbo  spin  echo  images.  -  Completed 

b)  Perform  semi-quantitative  and  quantitative  analysis  of  ex  vivo  prostate  samples 

We  reported  on  an  interim  statistical  analysis  during  our  last  annual  report.  We  just  finished  accruing  data 
from  our  final  ex  vivo  patients  this  past  month  and  are  in  the  process  of  compiling  a  final  statistical 
assessment  of  our  ex  vivo  data.  This  will  be  completed  by  the  end  of  this  current  NCE  period. 

c)  Statistically  analyze  MR-based  images  and  pathological  metrics 

Data  generation  continues  to  be  in  process  (MR  images  &  Pathology  metrics).  Figure  4  shows  a  typical 
panel  of  images  we  create  for  each  prostate  imaged.  We  also  have  the  ability  to  produce  a  variety  of  MR- 
EPT  images  based  on  different  algorithms  we  developed  (and  described  in  previous  reports)  (Fig  5).  We 
have  fully  reconstructed  47  of  our  50  cases. 
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Figure  4:  Typical  panel  of  data  images  recorded  and  to  be  used  for  statistical  analysis,  o  denotes  conductivity  images. 
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Figure  5.  Typical  panel  of  conductivity  images  based  on  applying  different  MR-EPT  algorithms  we  have  developed.  This  particular  example 
displays  a  single  slice  from  a  patient.  Five  different  reconstruction  algorithms  (discussed  in  previous  reports)  were  used  to  produce  these 
conductivity  images.  We  plan  to  statistically  analyze  each  of  the  algorithms  when  evaluating  the  clinical  efficacy  of  MR-EPT. 


Task  2  Milestones: 

1.  Assessment  of  the  clinical  potential  MR-EPT  combined  with  multi-parametric  MR  might  have  for 
prostate  imaging  -  in  process  (TBC) 

2.  1  peer-reviewed  publication  submitted  -  in  process  (TBC) 


SPECIFIC  AIM  3:  TO  EVALUA  TE  MR-EPT  IN  AN  IN  VIVO  COHORT  OF  PA  TIENTS 
Major  Task  1:  Evaluate  in  vivo  MR-EPT  and  multi-parametric  MR  imaging 

a)  Record  MR-EPT  and  multi-parametric  MR  sequences  of  in  vivo  and  ex  vivo  prostates 

Recruiting  patients  for  our  in  vivo  protocol  is  in  progress.  We  have  enrolled  5  men  into  our  in  vivo  protocol. 
We  have  conducted  an  initial  analysis  of  our  in  vivo  data  and  have  noted  that  there  are  significant  image 
artifacts  that  are  present.  This  may  be  due  to  a  phase  wrapping  issue  or  a  low  signal  issue  associated  with 
the  much  larger  field  of  view  that  we  are  trying  to  image.  The  figure  on  the  right  shows  an  example  series 
of  phase  images  acquired  of  a  man’s  prostate.  The  individual  images  correspond  to  an  image  stack.  The 
unsmooth  variation  in  the  images  are  not  expected  nor  observed  in  ex  vivo  and  small  phantom  imaging. 
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Figure  6.  In  vivo  series  of  prostate  phase  images  acquired  from  a  volunteer.  Note  the  significant  phase-wrapping  artifact  present  in 
these  images.  This  artifact  produced  MR-EPT  conductivity  images  with  severe  artifacts. 


To  determine  is  this  is  due  to  the  large  field  of  view  or  to  the  tissue  being  sampled  in  vivo  we  recorded  data 
from  a  volunteers  calf.  Note  the  phase  images  on  the  right  do  not  show  the  significant  artifact  present  in 
the  in  vivo  phase  images  of  the  prostate.  Further,  the  conductivity  images  derived  from  these  phase 
images  clearly  show  the  outer  adipose  layer  and  the  tibia  and  fibula  bones  within  the  calf  musculature.  We 
are  in  the  process  of  evaluating  large  saline  tanks  to  determine  how  the  field-of-view  influences  the 
acquired  phase  images. 


Figure  7.  In  vivo  series  of  calf  images.  Note  that  the  phase  images  do  not  suffer  from  the  same  phase  wrapping  artifact  present  within 
the  in  vivo  prostate  images.  The  conductivity  images  (center  panel)  were  successfully  reconstructed  from  these  phase  mages  and  clear 
show  the  low  conductivity  tibia  and  fibula. 


We  have  conducted  larger  phantom  tests  with  this  MR  bore  to  better  simulate  in  vivo  conditions. 
Specifically,  we  created  gel  phantoms  in  a  5-gallon  bucket.  The  phantom  ended  up  consisting  of  a  gelatin 
cylinder  ~12”  in  diameter  and  12"  in  height.  A  gelatin  prostate  with  copper  sulfate  (MR  contrast)  and 
additional  salt  (to  provide  conductivity  contrast)  was  formed.  A  plastic  (low  conductivity)  bead  was  placed 
within  the  prostate  to  emulate  a  tumor.  Note  that  this  study  was  carried  out  specifically  to  try  to  explore  the 
phase  artifacts  we  observed  in  the  in  vivo  cases.  We  showed  with  this  phantom  that  the  phase  images  are 
artifact  free  if  we  set  the  appropriate  field  of  view. 
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Figure  8.  Large  gel  phantom  used  to  better  simulate  in  vivo  images  conditions.  A  large  cylindrical  gel  phantom  was  created;  a  smaller 
anatomically  accurate  gel  prostate  with  a  bead  inserted  (to  simulate  a  tumor)  was  embedded  within  the  large  cylindrical  gel  phantom  (left 
panels).  The  phase  images  acquired  show  a  very  uniform  phase  distribution  with  no  phase  wrapping  (central  panels).  The  reconstructed 
MR-EPT  conductivity  images  clearly  show  the  prostate,  tumor,  and  a  small  air  bubble  that  was  trapped  in  the  phantom  during  fabrication. 
This  experiment  was  used  to  explore  several  field  of  views  (FOVs)  and  imaging  directions  within  the  bore.  The  optimal  FOV  used  a  large 
region  in  which  the  coil  was  positioned  Anterior  to  Posterior  and  the  phase  direction  was  oriented  from  left  to  right. 


11 


Based  on  these  results  we  conducted 
addition  in  vivo  tests  on  a  normal  volunteer  to 
see  if  changing  the  FOV  helped  to  minimize 
phase  artifacts  and  enable  MR-EPT 
reconstruction  of  conductivity  images.  We 
explored  a  number  of  different  FOVs  and 
imaging  directions  (see  Table).  Figure  10  shows 
example  magnitude  and  phase  images  acquired  of  a 
volunteer  using  the  FOVs  and  directions  specified  in 
the  above  table.  Note  that  the  phase  images  do  not 
have  the  wrapping  artifact  that  was  present  in  our 
earlier  phase  images  of  in  vivo  prostate  (see  above). 
We  have  reconstructed  these  images  to  produce  the 
first  ever  in  vivo  MR-EPT  images  of  prostate  (Figure 
11). 

Direct  MR-EPT:  Average  Laplace 
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Figure  10.  In  vivo  series  of  MR-EPT  prostate  conductivity  images 
acquired  from  a  volunteer.  These  represent  the  first  ever  in  vivo 
conductivity  images  computed  using  MR-EPT  of  the  prostate.  Two 
different  reconstruction  approaches  are  presented  here:  the  direct 
average  Laplace  approach  and  our  inverse  approach  using  our 
newly  implemented  overlapping  strategy  to  smooth  internal 
boundaries.  Both  panels  of  images  clearly  show  the  prostate  with  a 
somewhat  heterogeneous  intraprostatic  region  and  a  low 
conductivity  peri-prostatic  region.  This  approach  will  be  used  for 
image  acquisition  and  reconstruction  in  our  remaining  in  vivo 
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Figure  9.  In  vivo  series  of  prostate  images  acquired  from  a 
volunteer.  These  panels  show  the  magnitude  (left)  and  phase 
(right)  images  acquired  using  the  different  FOVs,  coil  and  phase 
directions  as  specified  in  the  above  table.  Note  that  the  phase 
wrapping  present  in  Series  11  was  easily  unwrapped  since  it  was 
somewhat  continuous  in  the  angular  direction.  Series  11  provided 
the  best  conductivity  images  (see  Fig  11). 


b)  Perform  post-prostatectomy  pathological  assessment  of  extracted  prostates 
To  be  completed  during  the  remaining  NCE  period  of  this  program 

c)  Statistically  analyze  MR-based  images  and  pathological  metrics 
To  be  completed  during  the  remaining  NCE  period  of  this  program 

Task  1  Milestones: 

1.  Comparison  between  in  vivo  and  ex  vivo  MR-EPT  -  TBC 

2.  Preliminary  clinical  statistics  defining  utility  of  MR-EPT  combined  with  other  MR-imaging  variants  -  TBC 

3.  Initial  parameter  threshold  values  for  use  in  detecting  and  staging  prostate  cancer  -  TBC 

4.  1  peer-reviewed  publication  submitted  -  TBC 


KEY  RESEARCH  ACCOMPLISHMENTS 

•  Assembled  a  database  of  MR  and  MR-EPT  images  from  50  ex  vivo  human  prostate  samples;  this  will 
ultimately  be  made  available  to  the  research  community  once  we  have  completed  statistical  analysis  of 
this  database. 

•  Developed  a  novel  overlapping  scheme  that  enables  us  to  produce  MR-EPT  conductivity  images  with 
no  internal  boundary  artifacts. 

•  Developed  MR-EPT  permittivity  imaging  using  an  inverse  approach  similar  to  what  we  developed  for 
use  in  conductivity  imaging. 

•  Produced  the  first  ever  MR-EPT-based  conductivity  images  of  in  vivo  human  prostate  that  clearly 
shows  the  prostate 

REPORTABLE  OUTCOMES 
Manuscripts 

Borsic  A,  Perreard  I,  Mahara  A,  Halter  RJ,  “An  Inverse  Problems  Approach  to  MR-EPT  Image  Reconstruction,” 
IEEE  Transactions  on  Medical  Imaging,  accepted  with  minor  modifications  June  2015.  (Appendix  1  - 
published  manuscript) 

Perreard  I,  Borsic  A,  Mahara  A,  Halter  RJ,  “Towards  Magnetic  Resonance  -  Electrical  Properties  Tomography 
(MR-EPT)  for  Prostate  Imaging,”  Medical  Physics,  to  be  submitted  September  2016  (draft  in  progress) 

CONCLUSION 

It  is  a  daily  challenge  for  clinicians  to  determine  whether  a  man  recently  diagnosed  with  prostate  cancer  has 
aggressive  disease  requiring  immediately  radical  therapy  or  indolent  disease  requiring  a  more  passive  watchful 
waiting  or  active  surveillance  approach.  This  program  is  focused  on  developing  Magnetic  Resonance  - 
Electrical  Property  Tomography  (MR-EPT)  specifically  for  prostate  imaging.  Over  the  past  year  we  have 
continued  to  optimize  our  MR-EPT  algorithms  for  estimating  the  prostate’s  electrical  conductivity  given 
magnetic  field  phase  and  magnitude  images  acquired  using  custom  MR  sequences.  We  have  conducted  a 
series  of  experiments  to  demonstrate  that  in  vivo  prostate  imaging  is  possible.  Additional  simulations  and 
phantom  experiments  have  been  conducted  to  explore  the  influence  of  field-of-view,  coil  position,  and  phase 
direction  have  on  our  MR-EPT  images  and  to  validate  that  we  can  image  in  vivo  objects  with  dimension  similar 
to  those  expected  in  human  imaging.  We  have  completed  enrollment  of  our  ex  vivo  study  (have  recorded  data 
from  50  ex  vivo  prostates)  and  are  currently  in  the  process  of  performing  a  final  statistical  assessment  of  this 
data.  Finally,  the  first  ever  in  vivo  prostate  conductivity  images  that  clearly  show  the  prostate  and  surrounding 
tissues  were  generated  using  MR-EPT.  Over  the  course  of  the  final  NCE  period,  we  will  be  focusing  primarily 
on  clinical  data  acquisition  (in  vivo  cohorts)  and  in  analyzing  images  acquired  to  assess  the  potential  of  using 
MR-EPT  (coupled  with  other  MR  imaging  variants)  to  distinguish  between  aggressive  and  indolent  prostate 
cancer. 
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An  Inverse  Problems  Approach  to  MR-EPT 

Image  Reconstruction 

A.  Borsic*,  I.  Perreard,  A.  Mahara,  and  R.  J.  Halter 


Abstract — Magnetic  Resonance — Electrical  Properties  Tomog¬ 
raphy  (MR-EPT)  is  an  imaging  modality  that  maps  the  spatial 
distribution  of  the  electrical  conductivity  and  permittivity  using 
standard  MRI  systems.  The  presence  of  a  body  within  the  scanner 
alters  the  RF  field,  and  by  mapping  these  alterations  it  is  possible 
to  recover  the  electrical  properties.  The  field  is  time-harmonic, 
and  can  be  described  by  the  Helmholtz  equation.  Approximations 
to  this  equation  have  been  previously  used  to  estimate  conductivity 
and  permittivity  in  terms  of  first  or  second  derivatives  of  RF  field 
data.  Using  these  same  approximations,  an  inverse  approach  to 
solving  the  MR-EPT  problem  is  presented  here  that  leverages 
a  forward  model  for  describing  the  magnitude  and  phase  of 
the  field  within  the  imaging  domain,  and  a  fitting  approach  for 
estimating  the  electrical  properties  distribution.  The  advantages 
of  this  approach  are  that  1)  differentiation  of  the  measured  data 
is  not  required,  thus  reducing  noise  sensitivity,  and  2)  different 
regularization  schemes  can  be  adopted,  depending  on  prior  knowl¬ 
edge  of  the  distribution  of  conductivity  or  permittivity,  leading  to 
improved  image  quality.  To  demonstrate  the  developed  approach, 
both  Quadratic  (QR)  and  Total  Variation  (TV)  regularization 
methods  were  implemented  and  evaluated  through  numerical  sim¬ 
ulation  and  experimentally  acquired  data.  The  proposed  inverse 
approach  to  MR-EPT  reconstruction  correctly  identifies  contrasts 
and  accurately  reconstructs  the  geometry  in  both  simulations 
and  experiments.  The  TV  regularized  scheme  reconstructs  sharp 
spatial  transitions,  which  are  difficult  to  reconstruct  with  other, 
more  traditional  approaches. 

Index  Terms — Conductivity,  electrical  properties  tomography, 
inverse  problem,  magnetic  resonance,  permittivity,  primal 
dual — interior  point  method,  quadratic  regularization,  recon¬ 
struction,  total  variation  regularization. 


I.  Introduction 


Magnetic  resonance-electrical  prop¬ 
erties  TOMOGRAPHY  (MR-EPT)  is  an  imaging 
modality  that  maps  the  spatial  distribution  of  the  electrical  con¬ 
ductivity  and  permittivity  using  standard  clinical  MR  systems. 
This  technique  exploits  Bi  field  perturbations  associated  with 
the  objects  present  within  the  bore  of  the  scanner.  These  field 


Manuscript  received  July  07,  2015;  revised  July  27,  2015;  accepted  August 
03,  2015.  Date  of  publication  August  20,  2015;  date  of  current  version  De¬ 
cember  29,  2015.  Asterisk  indicates  corresponding  author. 

*A.  Borsic  is  with  NE  Scientific  LLC,  Lebanon,  NH  03766  USA  (e-mail: 
aborsic@ne-scientific .  com) . 

I.  Perreard  is  with  the  Department  of  Radiology,  Dartmouth  Hitchcock  Med¬ 
ical  Center,  Dartmouth  College,  Lebanon,  NH  03766  USA. 

A.  Mahara  and  R.  J.  Halter  are  with  the  Thayer  School  of  Engineering,  Dart¬ 
mouth  College,  Hanover,  NH  03755  USA. 

Color  versions  of  one  or  more  of  the  figures  in  this  paper  are  available  online 
at  http://ieeexplore.ieee.org. 

Digital  Object  Identifier  10. 1109/TMI.20 15. 2466082 


perturbations  can  be  used  to  estimate  the  spatial  distribution 
of  electrical  properties  within  the  objects  giving  rise  to  the 
perturbations. 

Early  work  in  MR-EPT  was  published  in  1991  [1],  but  the 
technique  has  only  recently  seen  a  resurgence  of  exploration 
due  to  the  clinical  potential  that  the  electrical  properties  may 
offer  in  terms  of  tissue  contrast.  Several  groups  have  focused  on 
developing  novel  image  reconstruction  methods  and  algorithms 
for  this  purpose. 

Reconstruction  algorithms  can  be  classified  as  direct 
methods ,  from  which  measured  B\  information  is  directly  used 
to  map  the  electrical  properties  (EPs),  and  as  inverse  methods , 
in  which  EPs  are  estimated  by  fitting  a  model  to  the  measured 
data.  The  initial  work  by  Haake  et  al.  [1]  is  based  on  using  the 
Helmholtz  equation  for  describing  the  time  harmonic  magnetic 
field  H  at  radio  frequencies  used  in  MRI  imaging: 


V2H  =  -LiL02kH  + 


^x(Vxff) 

k 


(i) 


where  k  =  e  —  j(a/u),  and  uj  is  the  angular  frequency  of 
the  H  field,  a  the  electrical  conductivity,  e  the  electrical  permit¬ 
tivity,  and  ii  the  magnetic  permeability.  Assuming  that  a  and 
e  are  piecewise  constant  or  slowly  varying  (i.e.  \7k  «  0),  the 
second  term  on  the  right  hand  side  of  (1)  can  be  neglected.  Con¬ 
sidering  only  the  MRI-measurable  positive  circularly  polarized 
component  H+  of  the  RF  transmit  field,  and  considering  k  as 
isotropic,  one  obtains: 


-1  V2fT+ 

/X£j2  H+ 


(2) 


which  expresses  the  EPs  as  a  function  of  the  fT+  field  and  its 
second  derivative.  One  method  for  measuring  H+  amplitude 
is  through  use  of  double — angle  mapping  techniques  [2],  H+ 
phase  is  assumed  to  be  half  of  the  spin — echo  phase,  an  assump¬ 
tion  valid  when  one  transmit  and  one  receive  coil  are  used,  and 
the  sensitivity  patterns  of  the  two  coils  have  a  similar  spatial 
distribution  but  a  reverse  polarity  [3].  Typically,  in  biological 
tissues,  /i  is  considered  equal  to  the  permeability  of  free  space 
Fo- 

Despite  the  simplifying  assumptions  (2),  this  approach  has 
successfully  been  used  to  fit  layered  models  [  1  ]  and  to  produce 
electrical  properties  images  in  post-mortem  animals  [4].  While 
this  approach  is  feasible,  the  2nd-order  differentiation  required 
for  computing  the  Laplacian  is  sensitive  to  noise,  and  therefore 
not  desirable.  In  2009,  Katscher  et  al  [3]  used  the  Gauss  the¬ 
orem  to  propose  an  alternate  formulation  that  decreases  the  dif¬ 
ferentiation  from  2nd  to  1  st  order.  With  this  transformation,  the 
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conductivity  and  permittivity  are  proportional  to  surface  inte¬ 
grals  of  the  gradient  of  H+  phase  and  amplitude,  over  the  sur¬ 
face  of  an  arbitrary  block  of  pixels  where  the  electrical  prop¬ 
erties  are  assumed  to  be  approximately  constant.  This  approach 
has  been  successfully  demonstrated  in  in-vitro  and  in  in-vivo  ex¬ 
periments  [3],  [5],  [6]  and  exhibits  decreased  noise  sensitivity 
compared  to  previously  considered  approaches. 

The  assumption  that  EPs  are  constant  or  slowly  varying  also 
does  not  broadly  apply  to  imaging  of  biological  tissues,  and  the 
errors  arising  from  violating  this  assumption  are  analyzed  in 
detail  by  Seo  et  al  [7]. 

An  approach  based  on  measuring  two  sets  of  77+  data  with 
multi-channel  MRI  systems  was  proposed  by  Zhang  et  al.  [8]. 
This  approach  reduces  artifacts  resulting  from  fast  spatial  vari¬ 
ations  of  EPs,  but  requires  measuring  components  of  77+  in  the 
(x,  y )  plane  of  the  scanner,  which  is  impractical  in  clinical  ap¬ 
plications,  where  only  the  z  component  along  the  direction  of 
the  static  magnetic  field  is  normally  available. 

Sodickson  et  al  [9]  derived  a  general  formulation  to  MR-EPT 
that  does  not  make  assumptions  on  the  distribution  of  EPs  and 
that  is  suitable  for  multi-channel  systems.  In  this  approach  con¬ 
ductivity,  permittivity,  phase,  phase  derivatives,  magnetization, 
and  magnetization  derivatives  are  assumed  to  be  unknowns. 
This  approach  has  produced  promising  numerical  results,  how¬ 
ever,  it  requires  computing  spatially-dependent  second  deriva¬ 
tives  and  may  therefore  be  sensitive  to  noise  in  experimental 
applications.  An  approach  based  on  formulating  the  dependence 
of  EPs  on  77+  via  the  convection — reaction  equation  has  been 
developed  by  Hafalir  et  al  [10].  This  approach  does  not  make 
any  restrictive  assumption  on  the  distribution  of  EPs.  However, 
this  method  also  requires  first  and  second  derivatives  of  77+  be 
computed  and  is  therefore  potentially  sensitive  to  noise;  despite 
this,  it  has  been  demonstrated  numerically  and  experimentally 
to  provide  results  that  are  superior  to  standard  methods  based 
on  (2). 

Inverse  approaches  to  MR-EPT  have  been  also  considered. 
An  algorithm  termed  CSI-EPT  (Contrast  Source  Inver- 
sion-EPT)  has  been  proposed  by  Balidemaj  et  al  [11].  This 
algorithm  does  not  make  any  assumption  regarding  the  distri¬ 
bution  of  EPs  and  is  based  on  an  inverse  formulation,  where 
the  EPs  are  fit  to  the  data  using  a  contrast  source  approach.  The 
method  has  been  shown  in  numerical  simulations  to  produce 
accurate  and  detailed  reconstructions  of  a  pelvis  model.  Total 
Variation  [12]  regularization  has  also  been  adopted  in  order 
to  reduce  this  method’s  sensitivity  to  noise.  To  the  best  of 
our  knowledge  this  method  has  not  yet  been  demonstrated 
on  experimental  data.  A  model — based  based  approach  to 
conductivity  reconstruction  which  incorporates  regularization 
techniques  has  also  been  proposed  by  Ropella  et  al  [13].  This 
approach  is  based  on  inverting  the  Laplacian  in  (3)  in  the 
Fourier  domain,  and  has  demonstrated  better  image  quality 
compared  to  traditional  direct  approaches  in  phantoms  and  in 
in-vivo  data. 

In  this  manuscript  we  present  a  novel  MR-EPT  reconstruc¬ 
tion  algorithm  resulting  from  further  development  and  enhance¬ 
ment  of  initial  work  developed  by  the  authors  [14],  [15].  The 
algorithm  is  based  on  an  inverse  approach  applied  to  (2).  Con¬ 
ductivity  and  permittivity  are  treated  separately  as  in  [4],  and 


taken  as  parameters  to  be  fitted.  A  forward  model  is  developed 
to  link  electrical  parameters  to  the  77+  data.  A  Jacobian  ma¬ 
trix,  based  on  the  forward  model,  is  defined  and  used  for  up¬ 
dating  the  model  parameters.  Regularization  is  used  to  stabilize 
the  inversion  for  parameter  estimation.  This  inverse  approach  to 
MR-EPT,  based  on  fitting  77+  data  with  a  model,  has  the  general 
advantage  compared  to  direct  methods  based  on  (2),  and  with 
respect  to  [9],  [10],  of  not  requiring  differentiation  of  measured 
77+  data,  and  is  therefore  less  sensitive  to  noise.  The  approach 
is  demonstrated  to  successfully  reconstruct  noisy  numerical  and 
experimental  data.  Because  this  approach  casts  MR-EPT  recon¬ 
struction  in  a  well  established  inverse  problem  framework,  two 
well  known  regularization  techniques,  Quadratic  Regularization 
and  Total  Variation  regularization  are  implemented  to  reduce 
the  effects  of  noise  and  to  stabilize  the  inversion.  Quadratic 
Regularization  leads  to  smoother  reconstructed  images,  while 
Total  Variation  regularization  is  able  to  produce  sharper  im¬ 
ages.  This  method,  being  based  on  the  approximate  relationship 
(2),  which  might  result  in  artifacts  at  the  interface  of  highly  en¬ 
terogenous  boundaries  as  discussed  in  [7].  However,  this  ap¬ 
proach  is  common  in  MR-EPT,  and  has  been  demonstrated  to 
produce  meaningful  images  [3],  [5],  [6],  [16],  [17]. 

Algorithms  based  on  (2)  have  been  applied  to  breast  cancer 
detection,  showing  potential  of  MR-EPT  as  a  diagnostic  tool. 
In  this  context  prior  structural  information  available  from  T1 
and/or  T2  weighted  MRI  imaging  has  been  used  to  enhance 
the  reconstructed  EPs  by  weighting  differently  variations  along 
the  normal  and  tangential  direction  with  respect  to  the  expected 
prior  features  [16],  [17].  Our  approach  also  enables  incorpo¬ 
rating  prior  structural  information,  so  that  preferential  directions 
of  change  of  EP  can  be  embedded  into  the  regularization  (e.g. 
[18],  [19]).  In  general  we  believe  that  the  approach  developed 
here  has  advantages  over  algorithms  based  directly  on  (2)  since 
it  does  not  require  differentiation  of  77+ .  This  may  make  it  more 
suitable  for  clinical  applications  where  noise  is  present.  It  is  dif¬ 
ficult  to  offer  comparative  remarks  with  respect  to  recent  inverse 
algorithms  [9]— [1 1],  as  these  comparisons  will  likely  require  ap¬ 
plication  of  the  different  algorithms  on  particular  tests  cases. 

Reconstruction  of  EPs  with  the  inverse  approach  developed 
her  entails  considering  the  EPs  of  every  pixel  in  the  image  as  an 
unknown  to  be  estimated.  For  three-dimensional  datasets  this 
can  result  in  an  excessive  computational  burden.  We  describe 
methods  for  splitting  the  imaging  domain  that  significantly  re¬ 
duce  this  burden. 

In  Section  II  we  introduce  our  reconstruction  approach  and 
develop  a  forward  model.  In  Section  III  we  describe  an  in¬ 
verse  formulation  for  reconstructing  EPs,  and  in  Section  IV  we 
discuss  an  implementation  using  Quadratic  Regularization  and 
one  using  Total  Variation  regularization.  Section  V  reports  nu¬ 
merical  experiments,  Section  VI  discusses  the  computational 
burden  of  the  methods  developed  and  offers  methods  for  re¬ 
ducing  it,  Section  VII  reports  image  reconstruction  results  from 
physical  experiments.  Finally  concluding  remarks  are  offered. 

II.  Forward  Model 

The  MR-EPT  reconstruction  approach  developed  here  is 
based  on  an  inverse  formulation,  in  which  a  model  is  fit  to  mea¬ 
sured  data.  In  this  context,  the  relationship  linking  the  77+  field 
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data  to  the  electrical  conductivity,  a,  and  permittivity,  e,  defines 
the  forward  model.  These  relationships  can  be  approximated  as 
[4]: 

a  =  —  V20(tf+)  (3) 

jlQU 


and 


1  V2|  H+ 

\H+\ 


(4) 


Equation  (3)  can  be  recognized  as  the  Poisson  equation: 


V20  =  fioua 


(5) 


where  for  simplicity  0  denotes  0(iT+).  Equation  (4)  takes  the 
form  of  the  Helmholtz  equation: 


where  a  is  a  scalar  Tikhonov  factor  controlling  the  amount  of 
regularization  and  4/ (a)  is  a  regularization  functional.  These 
functionals  are  often  quadratic  and  involve  first  or  second 
differential  operators  [21],  but  many  different  forms  beyond 
simple  differential  operators  have  been  proposed  in  literature. 
Two  different  functionals  often  used  in  inverse  problems, 
Quadratic  Regularization  and  Total  Variation  Regularization, 
are  described  and  implemented  here. 

IV.  Implementation 

In  order  to  implement  an  MR-EPT  reconstruction  based  on 
(9),  a  regularization  functional  needs  to  be  defined,  a  numeric 
method  for  computing  0(a)  as  defined  by  (5)  and  (7)  must  be 
implemented,  and  the  Jacobian  matrix  of  the  mapping  a 
0(a)  needs  to  be  computed. 


V2A  +  n0cj2eA  =  0  (6) 

where  for  simplicity  A  denotes  \H+\. 

This  manuscript  focuses  on  reconstructing  conductivity  from 
(5).  Permittivity  can  be  reconstructed  from  (6)  with  similar 
methods,  which  are  not  developed  in  this  manuscript.  Boundary 
conditions  appropriate  to  the  physics  describing  the  forward 
model  need  to  be  specified  in  order  to  solve  (5).  This  approach 
ultimately  aims  to  match  a  measured  phase  0meas  with  a 
simulated  0  computed  from  (5);  Dirichlet  conditions  defining 
the  phase  on  the  domain  boundary  are  therefore  ideal  for  this 
case.  This  conditions  is  specified  as 

4>{r)  =  4>{r)meas  Vr  G  dV.  (7) 


A.  Quadratic  Regularization 

In  a  first  implementation  example,  a  classical  quadratic  func¬ 
tional  can  be  used  to  transform  (9)  into 

(Tree  =  argmin  || 0(<r)  -  4>meas\ |2  +  a\\La\\2  (10) 

where  a  is  discretized  on  the  same  pixel  grid  as  that  used  by  the 
MR  scanner  to  map  0meas  •  This  converts  a  from  a  continuous 
function  into  a  finite  vector  of  discretized  values.  The  matrix 
L  is  a  regularization  matrix,  which  we  have  chosen  to  be  the 
Laplacian  of  the  conductivity  distribution,  a  relatively  common 
choice  [22],  [23]. 

By  applying  the  Newton-Raphson  method  to  (10),  an  update 
equation  for  the  conductivity  can  be  derived  as 


where  r  is  a  point  in  space  and  dfl  is  the  boundary  of  the  imaging 
domain.  0meas  is  perfectly  matched  at  the  boundary  through 
use  of  (7),  and  it  will  be  matched  point-by-point  inside  the  do¬ 
main  by  appropriate  adjustment  of  the  a  distribution.  It  is  worth 
noting  that  for  any  value  of  measured  data,  0meas ,  the  condi¬ 
tion  (7)  is  compatible  with  (5)  [20],  and  therefore  a  solution  to 
(5)  always  exists  and  is  unique  for  those  boundary  conditions 
[20].  Equations  (5)  and  (7)  define  a  forward  model  for  MR-EPT 
conductivity  reconstruction  via  the  inverse  approach  developed 
here. 

III.  Inverse  Formulation 

The  forward  model  linking  0  to  a  can  be  used  to  establish  an 
estimate  of  the  a  distribution  by  fitting  the  phase  0  predicted 
by  the  model  to  the  measured  phase  0meas  •  As  an  example,  this 
fitting  can  be  optimized  in  the  least  squares  sense  as 

CTree  =  argmin \\(f)(cr)  -  4>meas\\2  (8) 

where  arec  is  the  spatial  distribution  of  the  electrical  conduc¬ 
tivity  reconstructed  by  fitting  0(a)  to  0meas.  This  fitting  is  ac¬ 
complished  by  adjusting  a  to  minimize  the  discrepancy  between 
0(a)  and  0meas  in  the  L2-norm  sense;  the  dependence  of  the 
model  predicted  phase  0  on  a  is  explicitly  shown  in  (8). 

To  account  for  noisy  phase  measurements,  a  regularization 
term  [2 1  ]  is  adopted  to  stabilize  the  inversion,  transforming  (8) 
into 

<rrec  =  argmin  \\<p(cr)  -  4>meas ||2  +  a^r(a)  (9) 


5a  =  ~[JTJ  +  aLTL]  1 

x  [JT  (4>(a)  -  <f>meas)  -  aLTL(a  -  a*)]  (11) 


where  J  is  the  Jacobian  matrix  of  the  mapping  a  i— 0(a),  and 
a*  is  an  initial  conductivity  distribution.  A  uniform  conductivity 
can  be  used  as  an  initial  starting  distribution  a*  and  updated 
using  (11),  as  arec  =  a*  +  Sa.  A  single  update  is  sufficient  in 
this  case  since  the  forward  model  is  linear  in  a.  As  a  result,  the 
Newton-Raphson  method  finds  the  solution  to  (10)  in  one  step. 
In  order  to  apply  (11)  one  has  to  compute  0(a)  and  the  Jacobian 
matrix  J. 

The  forward  problem  a  4  0(a)  in  three  dimensions  consists 
of  solving 


<920  <920  <920 

dx 2  ^  dy 2  dz 2  ^ 


(12) 


where  the  (x,y,x)  are  the  coordinates  in  the  axes  (x,  y,  z\ 
which  are  defined  to  be  aligned  to  the  main  axes  of  the  image 
stack.  The  Partial  Differential  (12)  can  be  easily  discretized  on 
the  image  grid  using  Finite  Difference  schemes  [24]  expressing 
the  partial  derivatives  as 


<920 

cj)(x  -\-  hx^j 

-  2 4>(x)  +  4>(x  -  hx ) 

(13a) 

dx 2 

hi 

920 

_  4>(y  +  hy) 

-  2 4>(y)  +  4>{y  -  hv) 

(13b) 

dy 2 

hi 

920 

_  <j>(z  +  hz) 

-  2 cj){z)  +  cf)(z  -  hz ) 

(13c) 

dz 2 

h2z 
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where  hx,hy,  hz  are  respectively  the  pixel  spacings  along  the 
x,  y ,  z  axes  and  the  points  (x  —  hx),  x ,  (, x+hx ),  ( y  —  hy ),  y,  (?/+ 
hy),  (z  —  hz),  z,  and  (2:  +  )  translate  into  the  indices  within  the 

vector  of  discrete  phases  (j).  As  is  standard  in  Finite  Differences, 
the  partial  differential  equation  (12)  can  be  transformed  into  a 
linear  system 


Acj)  =  b 


(14) 


where  A  is  a  matrix  derived  from  applying  (13)  to  the  different 
pixels  in  the  image,  b  is  a  right  hand  side  (RHS)  derived  from 
the  evaluation  of  the  term  yuoa  in  (12)  at  each  pixel  location  in 
the  image,  and  0  is  the  vector  of  computed  phase  values. 

Appropriate  boundary  conditions  need  to  be  applied  to  the 
linear  system  in  (14).  The  following  approach  is  used  here 

•  Boundary  Pixels:  Dirichlet  boundary  conditions  are  estab¬ 
lished  for  all  the  pixels  on  the  boundary  surface  of  the 
image  stack  by  placing  a  1  on  the  diagonal  element  A(i,i), 
where  i  is  the  index  of  each  of  these  pixels,  and  setting 

—  0(^)meas* 

•  Internal  Pixels:  for  all  the  pixels  internal  to  the  domain 
(defined  as  being  at  least  1  pixel  away  from  the  boundary), 
(13)  is  applied  and  the  RHS  is  computed  as  b(i)  =  yuja(i), 
where  i  is  the  index  of  the  considered  pixel. 

Solving  the  linear  system  (14)  results  in  the  vector  0  of  com¬ 
puted  phase  values. 

A  similar  approach  is  used  for  building  the  matrix  L ,  which 
is  a  Laplacian  operator.  The  purpose  of  L  is  to  limit  high  fre¬ 
quency  spatial  variations  in  the  image,  which  typically  arise 
from  noise  in  the  data.  The  term  ||Lcr||2  takes  large  values  for 
large  spatial  variations  or  for  high  frequency  spatial  variations 
in  cr,  which  penalizes  their  presence  in  the  reconstructed  image 
(10).  Elements  of  L  corresponding  to  the  interior  domain  are 
computed  using  the  same  finite  differences  applied  to  solve  the 
forward  problem  (13).  Elements  of  L  corresponding  to  the  do¬ 
main  boundary  are  defined  using  mirroring  boundary  condi¬ 
tions.  Specifically,  if  the  term  {x  +  hx)  in  (13a)  falls  outside  the 
domain,  (x  +  hx)  is  assumed  equal  to  (x  —  hx )  and  the  second 
derivative  for  that  image  location  is  expressed  as: 

d24>  _  -2<f>(x)  +  2<t>{x  -  hx) 

dx2  ~  h2x  '  1  ’ 


Similar  conditions  are  used  for  derivatives  in  y  and  z. 

The  Jacobian  matrix  in  (1 1)  can  be  computed  from  (14)  using 
the  following  matrix  identities 


90 

da 


djA^b)  _  ,  db 

da  da 


(16) 


where  £  is  the  derivative  of  b  with  respect  to  a.  In  this  case, 
£(i)  =  0  for  all  indices  i  corresponding  to  boundary  pixels  (for 
which  the  Dirichlet  condition  has  been  set)  and  £(z)  =  yuo  for 
all  the  indices  corresponding  to  interior  pixels  where  b(i)  is  set 
to  b(i)  =  yuja{i).  Equations  (14)  and  (16)  enable  calculation  of 
0  and  J,  which  are  in  turn  used  in  (1 1)  to  compute  arec. 

As  shown  later,  the  above  procedure  results  in  successful  re¬ 
constructions  on  synthetic  and  experimental  data.  The  use  of  the 
quadratic  regularization  functional  produces  conductivity  pro¬ 
files  that  are  relatively  smooth.  The  benefit  of  the  inverse  formu¬ 
lation  developed  here  (9)  is  that  different  functional  terms  can 


be  used  for  regularization.  In  the  next  subsection,  Total  Varia¬ 
tion  is  introduced  as  an  alternative  regularization  term. 


B.  Total  Variation  Regularization 

In  the  previous  sections,  a  framework  for  MR-EPT  image 
reconstruction  based  on  inverse  problems  was  developed.  One 
benefit  of  this  framework  is  that  different  regularization  terms 
can  be  chosen.  Regularization  functionals  affect  how  the  recon¬ 
structed  image  is  smoothed  and  different  choices  are  appropriate 
for  different  situations.  Total  Variation  (TV)  is  a  relatively  novel 
form  of  regularization  that  results  in  images  with  sharper  con¬ 
ductivity  transitions  as  compared  to  Quadratic  Regularization 
approaches  (e.g.  11).  Reconstructing  sharp  image  transitions 
can  be  challenging  in  MR-EPT.  In  direct  approaches  derivatives 
are  estimated  on  multiple  points  (e.g.  5,7,  or  9)  to  reduce  noise 
sensitivity,  but  this  results  in  smoothing.  In  the  inverse  approach 
described  above  the  regularization  smooths  the  reconstructed 
image,  again  for  reducing  sensitivity  to  noise.  The  use  of  TV 
regularization  instead  allows  reconstructing  fast  spatial  varia¬ 
tions  more  accurately. 

The  TV-based  inverse  formulation  for  (9)  is  expressed  as 


<jrec  =  argmin 


0  -  0r7 


TV{c 


(17) 


where  the  TV  functional  is  defined  as  TV  (a)  =  fQ  |V(cr)|dD 
and  Q  is  the  imaging  domain.  The  sharper  reconstructions  pos¬ 
sible  with  TV  regularization  arise  because  the  TV  functional 
remains  finite  for  step  changes,  while  quadratic  functionals  like 
fQ  |  V(cr) |2d£7,  or,  fQ  |V2(cr)|2dQ  (common  quadratic  regular¬ 
ization  functionals)  tend  towards  infinity.  As  a  result  of  the  large 
quadratic  functional  values  associated  with  fast  spatial  changes 
in  conductivity,  these  profiles  (i.e.  step  changes  in  conductivity) 
are  penalized,  rendering  the  reconstructions  smoother.  More  de¬ 
tailed  discussion  of  the  TV  function  properties  are  presented  in 
[12]. 

While  the  use  of  TV  is  desirable  for  reconstructing  sharp  vari¬ 
ations,  the  image  reconstruction  expressed  by  (17)  is  a  non-dif- 
ferentiable  optimization  problem,  and  special  techniques  need 
to  be  employed  to  minimize  1 1 0  (cr)  —  0meas  1 1 2  +  a  TV  (cr)  when 
acting  on  a.  In  this  case,  a  Primal  Dual-Interior  Point  framework 
developed  in  [25],  [26]  for  optimizing  (17)  was  used.  Specif¬ 
ically,  the  algorithm  named  “PD-IPM-L2-L1  Norm”  reported 
in  [25]  is  used.  The  MR-EPT  forward  model  0(cr)  and  Jacobian 
matrix  J  developed  in  Section  IV- A  are  input  into  the  optimiza¬ 
tion  algorithm,  resulting  in  arec  as  defined  in  (17). 

In  Sections  V  and  VII,  numerical  and  experimental  data 
are  used  respectively  to  compare  these  two  regularization 
approaches  and  demonstrate  how  the  smoothing/sharpening 
characteristics  of  the  reconstruction  can  be  tuned  by  the  regu¬ 
larization  functional. 


V.  Numerical  Experiments 

A  number  of  numerical  experiments  were  conducted  to  val¬ 
idate  this  inverse  formulation  for  MR-EPT  conductivity  recon¬ 
struction.  The  simulated  volume  of  interest  consists  of  a  cubic 
block  of  20  x  20  x  20  millimeters  split  in  half  by  passing  a 
plane  through  the  center  of  the  cube.  One  half  of  the  cube  is  set 
to  have  a  conductivity  of  1  Sm-1  while  the  other  half  is  set  to 
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Fig.  1.  Simulated  domain  used  for  the  numerical  experiment:  a  20  x  20  X 
20  mm  cube  with  a  conductivity  of  1  Sm- 1  on  one  side  and  of  2  Sm- 1  on  the 
other  was  generated  in  MATLAB.  The  simulation  of  the  MR  phase  within  the 
cube  was  computed  for  a  3  T  magnet  using  a  Finite  Difference  discretization 
with  1  mm  resolution,  using  the  modeling  approach  developed  in  Section  II. 

2  Sm_1,  as  shown  in  Fig.  1.  While  these  simulated  conductivi¬ 
ties  are  generally  larger  than  typical  physiological  values,  they 
serve  as  a  good  test  do  demonstrate  that  the  algorithms  can  re¬ 
construct  a  unit  step  change  from  1  to  2  Sm-1  at  the  interface 
between  the  two  halves  of  the  numerical  phantom. 

The  simulated  cube  was  discretized  into  lxlxl  mm 
volume  elements,  and  Finite  Differences  were  used  to  compute 
the  MRI  phase  using  equations  (12)  to  (14),  and  assuming  a  3 
Tesla  static  magnet  strength.  A  Dirichlet  boundary  condition 
of  4>{r)  =  0  Vr  G  <9D  was  assumed  on  the  boundary.  As  dis¬ 
cussed  in  Section  II,  Dirichlet  boundary  conditions  can  be  used 
for  matching  a  measured  phase  at  the  domain  boundary  while 
computing  the  forward  solution.  In  the  absence  of  boundary 
data,  the  condition  0(r)  =  0  Vr  G  dfl  is  a  practical  condition 
that  is  compatible  with  the  Poisson  equation  and  results  in  a 
unique  solution.  This  choice  does  not  alter  the  reconstructed 
conductivity,  which  depends  on  the  Laplacian  of  the  phase-the 
value  of  which  is  enforced  within  the  domain  by  the  Poisson 
equation  itself. 

Fig.  2(a)  shows  a  2D  cross  section  of  the  computed  synthetic 
phase;  Fig.  2(b)  shows  the  same  synthetic  phase  with  20%  addi¬ 
tive  Gaussian  noise,  as  discussed  later  and  as  used  in  the  recon¬ 
structions.  The  curvature  (i.e.  Laplacian)  of  the  phase  is  more 
pronounced  in  the  left  part  of  the  domain  due  to  the  higher  con¬ 
ductivity  (2  Sm-1)  within  this  region. 

The  simulated  phase  data  was  input  into  three  different  re¬ 
construction  algorithms.  The  first  is  an  implementation  of  the 
direct  approach  developed  by  Katscher  et  al  [3].  This  algorithm 
is  based  on  numerically  solving  (2)  and  using  a  volume  of  in¬ 
tegration  to  improve  robustness  to  noise.  Specifically,  the  order 
of  differentiation  is  decreased  from  2nd  to  1  st  order  using  the 
Gauss  theorem  to  convert  the  volume  integral  of  (2)  to  a  sur¬ 
face  integral  which  has  only  first  derivatives  of  the  field  vari¬ 
ables.  This  represents  an  approach  that  directly  computes  the 
output  conductivity  as  a  function  of  the  input  phase  distribution 
by  taking  1  st  order  derivatives  and  computing  their  integral  over 
specific  integration  volumes.  In  our  implementation  of  this  al¬ 
gorithm,  an  integration  volume  of  5  x  5  x  5  pixels  was  defined 
and  phase  derivatives  were  estimated  using  Savitzky-Golay  [27] 
filters  involving  7  points  (three  points  per  side  of  the  considered 
pixel).  Katscher  et  al.  uses  similar  integration  volumes,  and  5  to 
9  points  for  derivative  estimation.  We  found  that,  for  our  numer¬ 
ical  data,  using  only  7  points  represents  a  good  compromise  be¬ 


tween  image  sharpness  and  noise  sensitivity,  while  larger  num¬ 
bers  of  points  lead  to  smoother  images  with  limited  ability  to 
identify  sharp  transitions.  The  second  and  third  algorithms  are 
implementations  of  the  inverse  formulation  approach  developed 
in  this  manuscript  with  Quadratic  Regularization  (10)  and  with 
TV  regularization  (17),  respectively. 

Reconstructions  for  the  three  different  algorithms  are  shown 
in  Fig.  3  using  a  fixed  gray-scale  for  all  figures  spanning  conduc¬ 
tivity  values  from  0.5  to  2  Sm_1 .  The  top  row  of  the  figure  rep¬ 
resents  reconstruction  with  our  own  implementation  of  the  algo¬ 
rithm  proposed  in  [3],  the  second  row  represents  reconstruction 
with  the  inverse  algorithm  using  Quadratic  Regularization,  and 
the  last  row  reconstructions  with  the  inverse  algorithm  using  TV 
regularization.  Different  levels  of  simulated  noise  were  added  to 
the  phase  data  for  each  column  of  the  reconstructions.  Specif¬ 
ically,  noise  levels  of  5%,  10%,  15%,  and  20%  were  added  to 
columns  1  to  4,  respectively.  Noise  was  generated  by  extracting 
values  from  a  Gaussian  random  distribution.  These  values  were 
scaled  to  generate  a  particular  percent  noise  level  (5%  to  20%) 
and  added  to  the  noiseless  simulated  phase  of  Fig.  2.  As  an  ex¬ 
ample,  for  a  1%  noise  level,  an  input  noise  image  v  is  scaled 
such  that  (std(u)/std(0))  =  0.01,  where  std  indicates  the  stan¬ 
dard  deviation.  Fig.  2(b)  shows  an  example  2D  cross  section  of 
the  noisy  simulated  phase  for  a  20%  noise  level. 

For  the  two  inverse  algorithms  an  optimal  value  of  the 
Tikhonov  factor  was  found  empirically,  and  used  for  each  of 
the  different  noise  levels.  A  value  of  1  x  10-5  was  used  for 
the  Quadratic  Regularization  reconstruction,  and  a  value  of 
5  x  10-6  for  the  TV  reconstruction.  The  effect  of  the  choice  of 
the  Tikhonov  value  is  discussed  later  in  this  Section. 

The  direct  algorithm  successfully  reconstructs  the  phase 
information,  showing  a  higher  conductivity  on  the  left  side 
of  the  domain,  and  a  vertical  transition  between  the  more 
conductive  and  less  conductive  regions.  The  gray,  3-pixel 
wide,  band  present  at  the  boundary  of  the  image  is  an  artifact 
of  using  derivative  filters,  which  cannot  operate  in  proximity 
of  the  boundary;  the  derivative  filters  used  here  require  three 
pixels  per  side  of  the  considered  pixel.  The  inverse  quadratic 
algorithm  successfully  reconstructs  the  conductivity  profile 
showing  a  more  conductive  and  a  less  conductive  region  on 
the  left  and  right  sides,  respectively.  This  algorithm  provides 
a  similarly  smooth  transition  at  the  conductivity  boundary  as 
compared  to  the  direct  algorithm,  but  visually  is  more  stable 
in  the  presence  of  noise.  Total  Variation-based  reconstructions 
exhibit  a  sharp  transition  in  the  conductivity  distribution,  re¬ 
sulting  in  a  better  overall  estimation  of  the  original  data.  This 
algorithm  is  also  the  least  sensitive  to  noise  for  this  dataset. 
TV  is  an  appropriate  image  prior  for  distributions  with  sharp 
transitions  like  the  one  used  for  these  tests,  and  therefore  likely 
to  produce  better  results  compared  to  other  methods.  All  the 
images  produced  by  the  inverse  approach  present  a  border  of 
one  pixel.  Pixels  on  the  boundary  are  not  estimated  as  they  can 
be  affected  by  the  boundary  condition  </>(r)  =  </>(r)meas- 

Fig.  4  shows  the  effect  of  the  Tikhonov  factor  on  the  re¬ 
constructed  images.  This  figure  shows  2D  cross-sections  of 
three — dimensional  reconstructions  with  the  Quadratic  Reg¬ 
ularization  inverse  algorithm  for  a  noise  level  of  10%  and 
for  different  values  of  the  Tikhonov  factor,  which  are,  from 
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Fig.  2.  2D  cross  section  in  the  horizontal  plane  of  the  3D  computed  MRI  phase  for  the  conductivity  distribution  shown  in  Fig.  1 .  Subfigure  (A)  shows  the  noiseless 
phase,  and  subfigure  (B)  the  noisiest  phase  (20%  noise  level)  fed  into  the  numerical  simulations.  As  expected  the  computed  phase  has  a  higher  curvature  corre¬ 
sponding  to  the  more  conductive  region  (left  of  the  picture),  as  the  Laplacian  takes  larger  values  for  high  conductivity  regions  compared  to  less  conducting  regions 
(right  of  the  picture).  The  units  used  for  plotting  the  phase  are  radians,  and  the  range  is-0.08  to  0.00  radians. 
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Fig.  3.  2D  cross-sections  of  three-dimensional  reconstructions  of  synthetic  phase  data  with  noise  levels  of  5%,  10%,  15%,  and  20%,  associated  with  columns  1  to 
4,  respectively.  Each  row  in  the  figure  represents  a  reconstruction  with  a  different  algorithm:  the  top  row  shows  (for  different  levels  of  noise)  reconstructions  with 
our  own  implementation  of  the  algorithm  proposed  in  [3].  We  use  an  integration  volume  of  5  x  5  x  5  pixels,  and  estimate  derivatives  on  3  pixels  per  side  of  the 
central  pixel,  this  results  in  a  band,  3 -pixels  wide,  around  the  image  that  cannot  be  reconstructed.  The  reconstruction  successfully  identifies  the  conductivity  dis¬ 
tribution,  showing  a  left-to-right  transition.  The  second  and  third  row  of  the  figure  represent  reconstructions  with  the  developed  inverse  approach,  using  Quadratic 
Regularization  (QR)  and  Total  Variation  (TV)  regularization,  respectively.  The  QR  algorithm  is  able  to  identify  the  left — to — right  change  in  conductivity  and 
to  describe  it  with  a  certain  degree  of  smoothness,  which  is  characteristic  of  this  type  of  regularization.  The  TV  algorithm  is  able  to  identify  and  describe  with 
pronounced  sharpness  the  left — to — right  change  in  conductivity.  In  this  specific  case  the  TV  algorithm  also  appears  particularly  robust  to  noise.  For  both  QR  and 
TV  algorithms  an  optimal  Tikhonov  factor  was  found  empirically,  and  maintained  constant  across  the  different  levels  of  noise.  Specifically  a  value  of  1  x  10~5 
and  of  5  x  10“  6  was  used  respectively  for  the  two  algorithms.  With  the  parameters  used,  the  inverse  based  algorithms  seem  to  fare  batter  in  the  presence  of  noise 
compared  to  the  direct  algorithm,  as  it  would  be  expected  from  the  fact  that  they  do  not  need  to  differentiate  input  phase  data.  All  figure  are  represented  on  a 
grayscale  ranging  from  0  to  2  Sm"1. 


left  to  right,  1  x  10"7,  1  x  1(T6,  1  x  1(T4,  and  1  x  1(T3. 
These  values  bracket  the  optimal  value  of  1  x  10“ 5  used  in 
the  reconstructions  of  Fig.  3.  The  two  left  reconstructions  are 
under — regularized  (i.e.  a  too  small  Tikhonov  value  results 
in  noisy  images).  The  two  right  reconstructions  are  instead 
over — regularized  (i.e  a  too  large  value  of  Tikhonov  factor 


results  in  reconstructions  that  are  not  as  sensitive  to  noise,  but  in 
which  the  spatial  transitions  have  been  smoothed  excessively). 
The  value  of  1  x  10~5  used  in  Fig.  4  therefore  represents  a 
good  compromise  between  sensitivity  to  noise  and  resolution. 

Conductivity  profiles  across  the  middle  of  the  domain  and 
passing  through  the  transition  point  were  extracted  from  the 
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Fig.  4.  Effect  of  the  Tikhonov  factor  on  reconstruction.  The  four  figures  show  2D  cross-sections  of  three-dimensional  reconstructions  of  synthetic  phase  data  with 
a  fixed  noise  level  of  10%  and  a  varying  Tikhonov  factor  of  a  value  ofl  x  10-7,  1  x  10~6,  1  x  10-4,  and  1  x  10~3,  from  left  to  right.  The  chosen  values 
bracket  the  optimal  value  of  1  x  10“ 5 ,  which  was  found  empirically  and  used  in  the  reconstructions  in  Fig.  3.  The  values  1  x  10“ 7  and  1  x  10“ 6  result  in  an 
under-regularized  image  which  is  sensitive  to  noise.  The  values  of  1  X  10-4,  1  x  10“ 3  result  in  images  that  are  over-regularized,  where  sensitivity  to  noise  has 
been  greatly  reduced,  but  spatial  transitions  have  been  overly  smoothed.  All  figures  are  represented  on  a  grayscale  ranging  from  0  to  2  Sm_1 . 
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Fig.  5.  Subfigure  A  shows  a  plot  of  the  conductivity  value  on  a  horizontal  line  crossing  the  cube  of  Fig.  1,  where  the  vertical  axis  represents  the  conductivity 
value  and  the  horizontal  axis  the  position  along  the  left-to-right  direction  inside  the  conductivity  cube.  Subfigures  B,  C,  and  D  represent  a  similar  plot  for  the 
reconstructed  conductivity  values,  for  a  noise  level  of  10%,  respectively  for  our  own  implementation  of  the  algorithm  proposed  in  [3],  for  the  inverse  reconstruction 
with  Quadratic  Regularization,  and  for  the  inverse  reconstruction  with  Total  Variation  regularization.  In  these  three  subfigures  the  bold  continuous  line  represents 
the  reconstructed  values,  and  the  thin  dotted  line  the  true  conductivity  value,  as  in  subfigure  A.  The  algorithm  in  subfigure  B  and  C  have  a  similar  performance  in 
terms  of  how  steeply  they  can  describe  the  sharp  transition,  while  the  algorithm  in  subfigure  D  (Total  Variation)  is  able  to  describe  the  step  conductivity  change 
much  more  accurately. 


three  reconstructions  with  the  three  different  algorithms  for 
the  10%  noise  level  case,  and  are  shown  in  Fig.  5.  The  orig¬ 
inal  conductivity  profile  (Fig.  5(a))  used  for  generating  the 
synthetic  phase  data  shows  the  true  transition  each  algorithm 
attempts  to  recover.  The  bold  lines  depict  the  reconstructed 
profile,  while  the  true  conductivity  profile  is  shown  as  a  dotted 
line  in  Figs.  5(a),  5(b),  and  5(c).  While  the  direct  approach 
(Fig.  5(a))  and  the  inverse  approach  with  Quadratic  Regu¬ 
larization  (Fig.  5(b))  are  substantially  equivalent  in  terms  of 
how  fast  they  describe  the  transition  from  high  conductivity 
to  low  conductivity,  the  reconstruction  with  TV  regularization 
(Fig.  5(c))  shows  a  much  steeper  transition  in  the  reconstructed 
profile,  and  is  therefore  a  more  accurate  reconstruction  of  the 
step  conductivity  profile.  The  inverse  formulation  developed 
here  enables  one  to  select  regularization  functionals  that  are 
appropriate  for  the  problem  at  hand.  Besides  TV  and  similar 


edge-preserving  techniques  it  is  possible  to  envision  using 
ad-hoc  functionals  that  incorporate  prior  structural  information 
extracted  from  other  anatomic  MR  images  [18],  [19]. 

VI.  COMPUTATIONAL  BURDEN  CONSIDERATIONS 

While  the  developed  approach  offers  flexibility  in  terms  of 
the  regularization  functionals  available  and  does  not  require  dif¬ 
ferentiating  the  noisy  input  images,  it  does  require  a  higher  com¬ 
putational  load  compared  to  methods  based  on  direct  differen¬ 
tiation.  This  computational  burden  can  be  measured  in  terms  of 
number  of  computations  and  amount  of  working  memory  re¬ 
quired.  We  will  show  that  the  computational  burden  can  be  sig¬ 
nificant,  for  images  with  a  moderate  number  of  pixels,  but  that 
it  is  possible  to  reduce  it  by  subdividing  the  image  domain  and 
by  reconstructing  smaller  subdomains  individually.  We  will  use 


BORSIC  et  al.\  AN  INVERSE  PROBLEMS  APPROACH  TO  MR-EPT  IMAGE  RECONSTRUCTION 


251 


PHANTOM  PHOTOGRAPH 


t  L 


GELATIN  STRIPES  LOCATIONS  (LIGHT  COLOR  IN  MRI  IMAGE) 

Spatial  Location  [Pixel  Number] 

MRI  MAGNITUDE  IMAGE 

Fig.  6.  Top:  photographic  representation  of  the  conductivity  phantom.  A  con¬ 
ductivity  phantom  was  generated  by  slicing  a  slab  of  gelatin  in  slices  of  20  mm, 
15  mm,  10  mm,  and  5  mm  of  thickness.  The  slices  were  positioned  inside  a 
polymer  housing  and  secured  to  it  by  slightly  heating  the  housing  and  letting 
the  gelatin  in  contact  with  the  housing  walls  slightly  melt  before  re-solidifica- 
tion.  The  phantom  was  then  filled  with  a  saline  solution  (not  shown).  The  gelatin 
had  a  conductivity  of  1.8  Sm-1  and  the  saline  solution  of  4.1  Sm-1.  Copper 
sulfate  (CuS04)  was  added  to  the  gelatin,  so  that  the  gelatin  slabs  would  show 
with  a  different  intensity  in  traditional  MR  amplitude  images.  Bottom:  an  MR 
amplitude  image  shows  the  phantom.  Gelatin  slices  appear  in  a  brighter  gray  due 
to  the  inclusion  of  Copper  Sulfate  (note  that  gelatin  slices  have  a  lower  conduc¬ 
tivity  value  compared  to  the  saline  solution  (1.8  Sm-1  versus  4.1  Sm-1),  and 
therefore  they  appear  darker  in  the  MR-EPT  reconstructions  of  Fig.  7). 


this  approach  for  reconstructing  the  experimental  MRI  data  pre¬ 
sented  in  Section  VII. 

Image  reconstruction  with  the  developed  inverse  formulation 
requires  solving  the  forward  problem  (5),  using  for  example  a 
finite-difference  scheme,  and  updating  the  conductivity  or  per¬ 
mittivity  values  with  a  Gauss  Newton  update  as  in  (11).  For 
moderate  to  relatively  large  imaging  volumes  solving  (5)  or 
(11)  presents  challenges.  For  example,  an  imaging  volume  of 
256  x  256  x  256  pixels  requires  the  forward  problem  to  be  solved 
for  2563,  or  more  than  16  million  unknowns.  Problems  of  this 
size  can  be  efficiently  solved  on  a  desktop  computer  in  a  few 
minutes  of  computation  time  using  specialized  algorithms,  such 
as  those  based  on  Algebraic  Multigrid  Methods  [28].  These  al¬ 
gorithms  are  not  readily  available  in  computing  environments 
like  MATLAB,  and  require  the  user  to  procure  specialized  li¬ 
braries.  In  addition  to  the  time  required  for  forward  solving, 
there  is  a  particular  challenge  posed  by  the  memory  utilization 
required  to  solve  (11).  The  number  of  rows  and  columns  in  the 
Jacobian  matrix  J  is  equal  to  the  number  of  unknowns  in  the 
problem.  In  the  particular  case  of  a  256  x  256  x  256  imaging 
domain,  the  Jacobian  would  be  a  2563  x  2563  matrix,  requiring 


memory  allocation  beyond  the  limits  of  any  personal  computer. 
The  approach  we  have  implemented  subdivides  the  imaging  do¬ 
main  into  several  smaller  subdomains  over  which  the  Jacobian 
is  formed;  this  significantly  reduces  the  computational  burden 
and  memory  requirements.  For  a  TV  x  N  x  N  imaging  domain, 
the  size  of  the  forward  problem  is  N3  and  the  memory  required 
for  storing  the  Jacobian  is  proportional  to  N3  *  N 3  =  N 6 ;  di¬ 
viding  the  domain  into  a  small  number  of  subdomains  immedi¬ 
ately  reduces  the  computational  burden. 

If  the  domain  is  subdivided,  for  example  into  cubic  blocks 
smaller  than  the  full  domain,  the  boundary  condition  (7)  can  be 
applied  at  the  interface  between  the  different  blocks.  The  for¬ 
ward  problem  can  then  be  solved  (5)  within  each  subdomain. 
Since  (7)  can  be  applied  everywhere  using  the  measured  phase 
values,  the  forward  problem  can  be  solved  on  a  subdomain  in¬ 
dependently  from  the  neighboring  subdomains,  and  image  re¬ 
construction  carried  out  on  each  subdomain  independently  from 
others. 

The  only  potential  dependence  between  neighboring  pixels 
is  introduced  at  the  interface  between  two  imaging  subdomains 
if  the  regularization  matrix  L  in  (11),  or  in  the  equivalent 
Total  Variation  formulation  (17),  correlates  neighbouring 
pixels  across  the  two  different  subdomains.  In  our  current 
implementation,  we  do  not  regularize  pixels  across  different 
imaging  subdomains;  this  enables  us  to  run  fully  independent 
reconstructions  on  portions  of  the  full  domain.  For  real  MRI 
data  (see  Section  VII),  we  split  the  imaging  domain  using 
this  sub-domain  technique  to  reduce  the  computational  burden 
associated  with  the  large  number  of  pixels  within  the  MR 
image  stack.  The  number  of  pixels  recorded  during  a  standard 
MR-EPT  scan  are  significantly  more  than  those  used  in  our 
numerical  simulation  (Section  V). 

To  substantiate  the  above  discussion,  we  report  the  computa¬ 
tional  time  and  the  required  amount  of  memory  for  storing  the 
Jacobian  matrix  for  different  subdomain  sizes  used  for  recon¬ 
structing  the  phantom  in  Fig.  9.  The  original  MRI  image  (ac¬ 
quisition  details  described  in  Section  VII)  has  a  size  of  144  x 
144  x  144  pixels.  Within  this  volume  a  region  of  interest  (ROI) 
of61x55x5  pixels  was  chosen  to  capture  the  curved  detail  of 
the  phantom  and  to  exclude  the  boundaries  of  the  plastic  con¬ 
tainer  used  to  house  it.  A  magnitude  image  of  the  selected  ROI 
is  shown  in  Fig.  9(b).  The  slices  of  the  ROI  were  split  in  1,2, 
or  3  parts  along  the  vertical  and  horizontal  directions,  giving 
raise  to  a  number  of  smaller  subdomains  for  image  reconstruc¬ 
tion.  Precisely  we  used  the  following  configurations:  2x1,2 
x  3,  3  x  2,  and  3x3,  where  the  first  number  indicates  how 
many  subdivisions  were  used  along  the  vertical  direction  and 
the  second  number  indicates  the  number  of  subdivisions  along 
the  horizontal  direction.  For  example,  the  configuration  3x3 
results  in  9  subdomains  in  total,  where  slices  are  subdivided  with 
a  3  by  3  grid  (see  Fig.  10.  All  the  information  (5  slices)  was  used 
in  the  third  dimension. 

Table  I  reports  timing  and  memory  usage  information.  All 
computations  were  performed  on  a  workstation  with  an  Intel 
Xeon  3503  CPU  running  at  2.40GHz,  with  8  GB  of  memory, 
using  the  Windows  7  Ultimate — 64  bit  operating  system; 
algorithms  were  implemented  and  run  in  the  MATLAB  envi¬ 
ronment.  An  important  difference  exist  between  the  algorithm 
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Fig.  7.  MR-EPT  Reconstructions  of  MR  data  acquired  on  the  phantom  of  Fig.  6.  Subfigure  A  shows  a  “traditional”  MR-EPT  reconstruction  with  our  implemen¬ 
tation  of  the  algorithm  proposed  in  [3].  A  border  of  3  pixels  is  present  all  around  the  boundary  of  the  image,  as  derivatives  are  estimated  using  three  pixel  values 
per  each  side  of  the  current  pixel,  therefore  reconstruction  can  occur  only  for  pixels  that  are  at  a  distance  of  3  pixels  from  the  domain  boundary.  The  algorithm 
reconstructs  correctly  the  alternating  conductivity  values,  corresponding  to  the  high  conducting  saline  solution,  and  less  conducting  gelatin  slices.  Gelatin  slices 
appear  in  fact  darker,  as  opposed  to  the  lighter  gray  in  which  they  appear  in  the  MR  magnitude  image  in  Fig.  6(b),  as  the  contrast  mechanism  in  MR-EPT  is  based 
on  the  electrical  properties,  while  the  MR  amplitude  image  is  sensitive  to  the  presence  of  Copper  Sulfate.  Subfigure  B  shows  a  MR-EPT  reconstruction  using 
the  developed  inverse  formulation  with  Quadratic  Regularization,  and  Subfigure  C  a  MR-EPT  reconstruction  using  the  developed  inverse  formulation  with  Total 
Variation  regularization.  All  figure  are  represented  on  a  grayscale  ranging  from  0  to  4  Sm-1 . 


implementing  Quadratic  Regularization  and  and  the  algo¬ 
rithm  implementing  Total  Variation  regularization  (TV):  the 
quadratic  update  equation  (11)  is  linear  with  respect  to  a  (the 
Jacobian  does  not  depend  on  a),  and  therefore  the  term  [JT  J 
+  aLT L]-1  only  needs  to  be  computed  once,  and  can  be  ap¬ 
plied  later  to  any  subdomain  of  the  image-a  fast  matrix — vector 
multiplication.  The  Total  Variation  regularization  algorithm 
instead  requires  an  iterative  cycle  on  each  subdomain  of  the 
image  (we  use  a  fixed  number  of  iterations,  10  cycles  per  sub- 
domain),  resulting  in  longer  reconstruction  times  compared  to 


Quadratic  Regularization  (QR).  For  both  algorithms  computing 
was  performed  in  double  precision,  timing  and  memory  usage 
information  is  therefor  relative  to  an  8-byte  representation  of 
floating  point  values.  The  coarser  subdivision,  2x1,  requires 
61  seconds  to  compute  with  QR  and  685  seconds  with  TV, 
using  162  Megabytes.  These  timing  and  memory  requirements 
reduce  respectively  to  1.5  and  86  seconds  for  the  smaller  3x3 
subdomain  configuration,  and  to  7  MB  of  memory  for  Jacobian 
storage.  While  the  QR  algorithm,  for  smaller  subdomains, 
presents  a  computational  burden  that  is  not  dissimilar  to  that 
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Fig.  8.  Phase  fitting  process.  This  figure  shows  how  the  MR-EPT  inverse  reconstruction  algorithm  is  fitting  the  phase  data.  The  plot  on  the  left  shows  results 
relative  to  the  QR  algorithm,  and  the  plot  on  the  left  result  relative  to  the  TV  algorithm.  The  dash-dotted  line  at  the  top  of  the  plots  represents  the  computed 
phase  value  across  a  numerical  phantom  of  Fig.  6  corresponding  to  a  uniform  initial  distribution  of  conductivity.  The  solid  bold  line  at  the  bottom  of  the  plots 
represents  the  true  measured  phase  across  the  phantom.  The  dips  and  peaks  of  the  measured  phase  represent  phase  changes  corresponding  to  the  alternating  high-low 
conductivity  values  encountered  across  the  gelatin  stripes  and  saline  volumes  that  constitute  the  phantom  (see  Fig.  6(a)).  The  dash-dotted  line  at  the  bottom  of  the 
plots  represents  the  fitted  phase,  resulting  from  the  inversion  procedure.  Despite  reconstructed  images  looking  quite  different  for  the  QR  and  TV  algorithms,  the 
fitted  phases  show  minor  differences.  This  is  to  be  expected,  as  different  values  of  reconstructed  conductivity  can  stem  from  minor  curvature  changes  in  the  phase. 
For  both  algorithms  the  fitted  phase  follows  closely  the  general  trends  of  the  measured  phase;  regularization  techniques  used  in  the  algorithm  help  to  ensure  that 
the  model  does  not  “overfif  ’  the  small  scale  phase  perturbations  that  represent  noise  in  the  data  (e.g.  small  dip  near  pixel  position  20). 
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Fig.  9.  Curved  phantom  photographic  image  (left)  and  MRI  magnitude  image 
(right).  This  phantom  was  build  with  a  similar  method  to  the  phantom  shown 
in  Fig.  6,  by  creating  a  gelatin  slab  which  is  immersed  in  a  saline  solution.  The 
gelatin  slab  has  a  conductivity  of  of  1.8  Sm-1  and  Copper  sulfate  (CuS04) 
was  added  to  it,  to  generate  contract  in  the  MRI  magnitude  image.  The  saline 
solution  has  a  conductivity  of  4.1  Sm-1 .  In  part  (A)  of  the  figure  is  indicated 
a  white  edge  that  approximates  the  perimeter  of  a  region  of  interest  which  was 
selected  for  image  reconstruction.  Part  (B)  of  the  figure  shows  the  MRI  magni¬ 
tude  corresponding  to  the  region  of  interest  used  in  the  reconstructions  of  Fig.  7. 

of  direct  approaches,  TV  reconstruction,  as  implemented  in  our 
algorithm,  still  presents  a  burden  that  is  significant.  Use  of  TV 
regularized  reconstruction  is  therefore  a  tradeoff  between  the 
benefits  offered  by  this  type  of  functional  and  the  time  con¬ 
sumed  in  the  reconstruction. 

VII.  Physical  Experiments 

The  proposed  MR-EPT  reconstruction  approach  is  demon¬ 
strated  on  two  phantom  datasets  acquired  using  a  Philips 
Achieva  3  T  MR  scanner  using  a  Philips  SENSE  Flex  S 
transmit/receive  coil.  A  first  conductivity  phantom  was  pre¬ 
pared  by  slicing  a  12  cm  x  7  cm  x  5  cm  block  of  gelatin  (a 
=  1.8  Sm-1)  into  slabs  of  different  thicknesses  (see  Fig.  6). 
Specifically,  slice  thickness  of  20  mm,  15  mm,  10  mm,  and  5 
mm  were  formed  and  placed  into  a  polymer  housing.  The  slices 


TABLE  I 

Computational  Time  For  Different  Subdomain  Configurations 


Config. 

Time  QR  Recon  [s] 

Time  TV  Recon  [s] 

Jac.  Mem.  [MB] 

2x1 

50.9 

684.8 

162 

2x2 

8.8 

232.9 

39 

3x2 

3.5 

133.3 

16 

3x3 

1.5 

85.9 

7 

were  positioned  inside  the  housing  and  adhered  to  the  walls  by 
slightly  heating  the  container;  the  heating  causes  a  thin  layer  of 
gelatin  to  melt  and  re-solidify  providing  adhesion.  Gelatin  was 
prepared  by  using  10%  by  weight  of  dry  porcine  gelatin,  and 
approximately  1%  by  weight  of  NaCl  to  adjust  conductivity  to 
the  desired  level.  A  saline  solution  with  a  conductivity  of  4.1 
Sm-1  was  produced  using  approximately  2%  NaCl  in  weight 
and  used  to  fill  the  gaps  between  the  slices  of  gelatin  (not  shown 
in  Fig.  6),  to  provide  an  inclusion  (gelatin)  to  background 
(saline)  contrast.  In  addition,  approximately  1%  by  weigh  of 
copper  sulfate  (CUSO4)  was  added  to  the  gelatin,  but  not  to  the 
saline,  in  order  to  provide  MR  contrast. 

Standard  MR  images  (e.g.  T1 -weighted  and  T2-weighted) 
were  not  expected  to  detect  the  electrical  properties  contrast  be¬ 
tween  the  gelatin  and  the  saline  solutions;  the  CuSCN  provides 
MR  contrast  so  that  the  structure  of  the  phantom  can  be  appre¬ 
ciated  in  standard  MR  images.  MRI  magnitude  images  depict  a 
high  intensity  (white)  where  the  gelatin  slabs  are  located  (corre¬ 
sponding  to  the  high  concentration  of  CuSCN),  while  the  saline 
solution  appears  as  a  low  intensity  (black)  region  (see  Fig.  6). 
This  magnitude  MR  image  provides  a  high  resolution  image  of 
the  experimental  configuration  to  which  reconstructed  MR-EPT 
images  can  be  compared.  A  standard  Spin  Echo  (SE)  MRI  se¬ 
quence  was  utilized  for  acquiring  the  magnitude  and  phase  data, 
with  the  following  settings:  resolution  of  80  x  80  x  20  pixels, 
field  of  view  (FOV)  of  160  x  160  x  60  mm,  repetition  time 
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Fig.  10.  2D  cross-sections  of  three-dimensional  reconstructions  of  the  phantom  of  Fig.  7.  This  figure  demonstrates  the  ability  of  the  inverse  algorithms  to  recon¬ 
struct  curved  boundaries,  and  also  serves  to  demonstrate  the  possibility  of  speeding  up  image  reconstruction  by  splitting  the  image  domain  in  smaller  subdomains 
to  be  reconstructed  independently.  The  left  column  represents  reconstructions  with  the  inverse  algorithm  with  Quadratic  Regularization  and  the  right  column 
reconstructions  using  Total  Variation  regularization.  Each  row  represents  a  different  subdomain  splitting  for  image  reconstruction,  as  discussed  in  Section  VI.  In 
the  top  row  MRI  slices  have  been  split  in  two  parts  vertically,  indicated  as  2  X  1  splitting,  in  the  second  row  the  slices  have  been  split  in  two  parts  vertically 
and  horizontally,  indicated  as  2  x  2  splitting,  and  the  third  and  fourth  represent  respectively  3x2  and  3x3  splittings.  Splitting  the  original  image  domain  in 
smaller  parts  allows  for  a  faster  reconstruction  as  discussed  in  Section  VI  and  reported  in  Table  I.  Some  discontinuities  are  present  at  those  image  locations  where 
two  different  subdomains  meet.  This  results  from  the  fact  that  different  subdomains  are  treated  as  separate  by  the  reconstruction  and  no  continuity  is  enforced.  In 
future  work  we  intend  to  introduce  some  correlation  between  neighbouring  subdomains,  through  regularization  techniques,  which  should  reduce  or  eliminate  these 
discontinuities.  In  both  sets  of  reconstructed  images,  though  some  minor  discontinuities  are  present  at  the  interface  between  subdomains  both  and  inverse  QR  and 
inverse  TV  algorithms  are  able  to  correctly  reproduce  the  curved  interface.  All  figures  are  represented  on  a  grayscale  ranging  from  0  to  4  Sm-1 . 


(TR)  —  600  ms,  and  echo  time  (TE)  —  7.79  ms.  Four  tem¬ 
poral  averages  were  used  in  order  to  improve  the  SNR  for  a  total 
acquisition  time  of  389  seconds.  MR  phase  information  was 
captured  together  with  amplitude  information  and  used  for  re¬ 
constructing  conductivity  images  with  the  inverse  implementa¬ 
tion  and  with  the  direct  method  proposed  by  Katscher  et  al.  [3]. 
For  the  inverse  formulation,  both  Quadratic  Regularization  and 
TV-based  regularization  approaches  were  used.  For  all  recon¬ 
structions,  data  extracted  from  a  central  portion  of  the  phantom, 


consisting  of  a  58  x  22  x  7  pixel  volume,  was  used.  All  three 
reconstruction  approaches  exhibit  regions  of  high  and  low  con¬ 
ductivity  corresponding  to  the  saline  and  gelatin,  respectively 
(see  Fig.  7).  The  same  grayscale  was  selected  for  each  of  the 
different  algorithms,  ranging  from  0  Sm-1  (black)  to  4  Sm-1 
(white).  It  is  important  to  note  that  gelatin  slices  present  a  lower 
conductivity  (1.8  Sm-1)  compared  to  the  saline  solution  (4.1 
Sm-1),  making  the  gelatin  slices  appear  in  a  darker  gray  level. 
This  is  in  contrast  to  MR  magnitude  images  (Fig.  6),  where 


BORSIC  et  al. :  AN  INVERSE  PROBLEMS  APPROACH  TO  MR-EPT  IMAGE  RECONSTRUCTION 


255 


gelatin  slices  appear  brighter  due  to  the  presence  of  (CUSO4). 
These  contrast  differences  observed  between  MR  magnitude 
images  (Fig.  6)  and  MR-EPT  reconstructions  (Fig.  7)  help  to 
validate  MR-EPT’ s  dependence  on  conductivity  and  not  on  typ¬ 
ical  MR  contrast  mechanisms  (e.g.  CUSO4). 

All  three  reconstruction  approaches  correctly  identify  the 
location  of  the  gelatin  strips  within  the  phantom  (Fig.  7).  The 
reconstruction  using  our  implementation  of  the  direct  method 
developed  by  Katscher  et  al.  [3]  has  a  border  of  3  pixels  over 
which  the  conductivity  values  are  not  reconstructed.  This 
artifact  stems  from  using  Savitzky-Golay  filters  with  3  pixels 
on  each  side  of  the  central  pixel  of  interest  to  estimate  the 
derivatives  at  these  locations.  In  the  inverse  formulation,  a 
single  pixel  around  the  border  is  used  to  match  the  measured 
phase  with  the  Dirichlet  boundary  condition  (7)  and  is  not 
reconstructed;  all  other  interior  pixels  are  fit  to  the  data.  In  this 
particular  experimental  configuration,  the  reconstruction  with 
Quadratic  Regularization  (Fig.  7(b))  provides  slightly  more 
contrast  as  compared  to  the  reconstruction  based  on  integration 
of  first  derivatives  (Fig.  7(a)).  Much  more  noticeable  is  the 
difference  between  the  reconstruction  based  on  Total  Variation 
(TV)  (Fig.  7(c))  regularization  and  the  other  two  reconstruction 
approaches.  TV  represents  a  good  prior  for  describing  the  sharp 
conductivity  transitions  that  are  present  at  the  gelatin-saline  in¬ 
terface;  using  this  form  of  regularization  significantly  enhances 
the  reconstructed  images.  The  reconstructed  conductivity 
(Fig.  7(c))  geometrically  aligns  with  the  phantom  configuration 
(as  viewed  visually  and  in  MR  magnitude  images)  in  that  the 
high  conductivity  gelatin  regions  (dark  stripes)  correspond 
to  high  intensity  MR  magnitude  regions  (light  stripes).  The 
TV-based  regularization  accurately  identifies  the  steep  transi¬ 
tions  that  are  characteristic  of  the  phantom,  while  the  direct 
and  quadratic-based  inverse  formulation  do  not  recover  the 
steep  transitions.  Both  quadratic  and  TV-regularized  inverse 
approaches  to  reconstruction  (Fig.  7(b)  and  (c))  exhibit  a  con¬ 
ductivity  discontinuity  through  the  center  of  the  image.  This 
discontinuity  arises  from  having  subdivided  the  reconstruction 
domain  into  an  upper  and  lower  subdomain  for  the  purpose 
of  speeding-up  reconstruction  as  discussed  in  Section  VI.  In 
the  current  implementation  of  this  approach  the  different  sub- 
domains  in  which  an  image  can  be  divided  are  reconstructed 
in  a  completely  independent  manner  (i.e.  we  do  not  introduce 
correlation  between  neighbouring  pixels  at  the  interface  of  the 
subdomains  in  the  regularization  matrices),  and  results  in  small 
discontinuities  observed  in  the  images.  Using  a  scheme  that 
incorporates  subdomain  overlap  as  presented  in  [29]  will  be 
explored  in  the  future  to  reduce  or  eliminate  this  artifact. 

The  fitting  process  resulting  from  the  inverse  approach  to  re¬ 
construction  can  be  evaluated  by  observing  how  the  fitted  phase 
progresses  from  the  initial  guess  to  the  final  reconstructed  phase 
(Fig.  8).  In  this  case,  the  estimated  phase  (y-axis)  is  plotted  as 
a  function  of  distance  (x-axis)  along  a  trajectory  passing  hori¬ 
zontally  through  the  striped  phantom  for  both  the  QR  and  TV 
algorithms.  The  dips  and  peaks  in  the  phase  data  are  caused 
by  the  varying  conductivity  values  across  the  phantom  gelatin 
slabs  and  saline  volumes.  True  measured  phase  is  indicated  with 
a  solid  line,  and  the  dash — dotted  line  at  the  top  of  the  plots 


represents  the  computed  phase  for  an  initial  uniform  conduc¬ 
tivity  distribution  of  1  Sm-1,  while  the  dash — dotted  line  at  the 
bottom  of  the  plots  represents  the  computed  phase  after  the  fit¬ 
ting  is  complete.  The  computed  phase  closely  tracks  the  mea¬ 
sured  phase,  with  the  exception  of  small-scale  features  associ¬ 
ated  with  noisy  measurements.  The  regularization  helps  to  en¬ 
force  an  accurate  fit  of  the  measured  phase,  but  to  avoid  fitting 
small  scale  variations  that  are  typically  associated  with  mea¬ 
surement  noise  and  errors. 

A  second  conductivity  phantom  was  prepared  using  the  same 
method  described  above  for  the  striped  phantom  of  Fig.  6.  In 
this  case  though  a  portion  of  a  round  slab  of  gelatin  was  cut  and 
placed  in  a  corner  of  a  polymer  housing  as  shown  in  Fig.  9.  This 
second  physical  phantom  was  used  to  evaluate  the  inverse  re¬ 
construction  algorithms  on  a  curvilinear  geometry.  The  striped 
phantom  of  Fig.  6  is  well  suited  for  reconstruction  with  TV 
regularization  algorithms,  as  it  presents  straight  step  changes, 
aligned  with  one  axis  of  the  image,  that  can  be  accurately  recon¬ 
structed  by  TV  algorithms.  The  curved  boundary  of  the  phantom 
in  Fig.  9  poses  a  challenge  to  TV  based  algorithms,  as  they 
are  know  to  tend  to  reconstruct  round  boundaries  with  a  stair¬ 
case  appearance.  The  reconstructions  of  Fig.  10  demonstrate 
that  both  the  QR  and  TV  algorithms  are  able  to  capture  and  re¬ 
produce  the  curved  nature  of  the  boundary  between  gelatin  and 
water.  Both  algorithms  exhibit  some  artifact  at  the  interface  be¬ 
tween  subdomains,  since  we  have  not  enforced  any  contiguity 
between  the  reconstructed  values  from  different  subdomains. 
TV  regularization  seems  in  this  case  to  suffer  less  from  the  stair¬ 
case  effect  compared  to  application  in  other  tomographic  appli¬ 
cations.  For  example  in  Electrical  Impedance  Tomography  the 
staircase  effect  seems  to  be  more  pronounced  [25].  We  believe 
this  might  be  due  to  the  fact  that  in  MR-EPT  data  is  measured 
everywhere  in  the  domain  and  not  only  at  the  boundary  as  in 
other  techniques.  This  might  help  to  drive  the  reconstruction  to¬ 
wards  more  realistic  results. 

VIII.  Conclusions 

A  novel  approach  to  MR-EPT  image  reconstruction  based 
on  an  inverse  formulation  has  been  developed.  The  approach  is 
based  on  the  observation  that  the  central  equations  of  MR-EPT 
can  be  inverted  and  used  to  describe  a  forward  model,  which 
in  turn  can  be  used  in  an  inverse,  data  fitting  approach  to 
MR-EPT  reconstruction.  This  approach  is  valid  for  reconstruc¬ 
tion  of  conductivity  from  measured  phase  information  and  for 
reconstruction  of  permittivity  from  measured  77+  amplitude 
information,  though  in  the  present  manuscript  we  develop 
and  demonstrate  this  approach  only  for  reconstruction  of  con¬ 
ductivity.  A  forward  model  for  computing  phase  information 
from  conductivity  data  is  presented  (Section  II).  In  addition,  an 
inverse  approach  based  on  Quadratic  Regularization  and  Jaco¬ 
bian  computation  has  been  developed  for  solving  the  MR-EPT 
image  reconstruction  problem  (Section  IV- A).  An  inverse  for¬ 
mulation  using  Total  Variation  as  a  functional  for  regularization 
has  been  developed  (Section  IV-B);  this  approach  enables  the 
reconstruction  of  sharper  conductivity  transitions  within  the 
imaging  domain.  Numerical  experiments  were  conducted  to 
validate  the  developed  inversion  approaches  in  the  presence 
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of  synthetic  noisy  data  (Section  V).  A  method  for  splitting 
the  imaging  domain  into  subdomains  has  been  developed  to 
reduce  the  computational  burden  arising  from  the  proposed 
inverse  approach  (Section  VI);  this  implementation  requires 
less  memory  utilization  and  fewer  computations  leading  to 
faster  reconstructions.  Subdividing  the  imaging  domain  into 
smaller  subdomains  results  in  significant  performance/memory 
utilization  gains,  as  both  the  memory  and  computation  time  are 
linked  to  the  size  of  the  Jacobian  matrix  (which  grows  with 
the  6th  power  of  the  size  of  the  imaging  domain).  Successful 
reconstructions  were  obtained  on  laboratory  data  collected 
from  a  phantom  experiment  that  included  alternating  regions 
of  high  and  low  conductivity,  and  from  a  phantom  with  a 
curved  boundary  (Section  VII).  Both  the  inverse  approach 
using  Quadratic  Regularization  and  the  approach  using  Total 
Variation  regularization  accurately  identified  the  position  of  the 
gelatin  slabs;  TV  regularization  more  accurately  reconstructed 
the  steep  conductivity  transitions  present  at  the  gelatin-saline 
interface  of  the  phantoms  imaged.  Benefits  of  this  reconstruc¬ 
tion  method  include:  1)  no  differentiation  is  required  on  the 
input  data,  therefore  decreasing  sensitivity  to  noise  compared 
to  other  methods  described  in  the  literature;  2)  different  reg¬ 
ularization  functionals  can  be  implemented,  depending  on  the 
expected  distribution  of  the  parameters  to  be  reconstructed. 
This  capability  is  particularly  useful  in  the  context  of  recon¬ 
structing  biomedical  data,  in  such  cases  custom  regularization 
functionals  can  be  constructed  ad-hoc  to  incorporate  prior 
anatomical  information  and  ultimately  enhance  the  quality  and 
robustness  of  the  reconstructed  images. 
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